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We propose a simple algorithm for generating normally distributed pseudo random numbers. 
The algorithm simulates N molecules that exchange energy among themselves following a simple 
stochastic rule. We prove that the system is ergodic, and that a Maxwell like distribution that 
may be used as a source of normally distributed random deviates follows in the N — > oo limit. 
The algorithm passes various performance tests, including Monte Carlo simulation of a finite 2D 
Ising model using Wolff's algorithm. It only requires four simple lines of computer code, and is 
approximately ten times faster than the Box-Muller algorithm. 
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Pseudo random number (PRN)generation is a subject 
of considerable current interest M. Deterministic algo- 
rithms lead to undesirable correlations, and some of them 
have been shown to give rise to erroneous results for ran- 
dom walk simulations 0, Monte Carlo (MC) calculations 
{§,§), and growth models ||. Most of the interest has 
been focused on PRN's with uniform distributions. Less 
attention has been paid to non-uniform PRN generation. 

Sequences of random numbers with Gaussian probabil- 
ity distribution functions (pdf 's) are needed to simulate 
on computers gaussian noise that is inherent to a wide 
variety of natural phenomena || . Their usefulness tran- 
scends physics. For instance, numerical simulations of 
economic systems that make use of so called geometric 
Brownian models (in which noise is multiplicative) also 
need a source of normally distributed PRN's 0. There 
are several algorithms available for PRN's with Gaus- 
sian pdf's §. Some, such as Box-Muller's algorithm, 
require an input of uniform PRN's, and their output of- 
ten suffers from the pitfalls of the latter Q . Robustness 
is therefore a relevant issue. In addition, Box-Muller's 
algorithm is slow and can consequently consume signif- 
icant fractions of computer simulation times p0| . The 
comparison method demands several uniform PRN's per 
normal PRN, and is therefore also slow ElJ]. Use of ta- 
bles [|l2| is not a very accurate method. Algorithms that 
are related, but not equivalent, to the one we propose 
here have been published (h],[l| , but they are somewhat 
cumbersome to use. In addition, no proof of their validity 
has been given. 

We propose here a new algorithm for the generation 
of normally distributed PRN's that is quite simple and 
fast. It is a stochastic caricature of a closed classical 
system of N particles. Their velocities provide a source 
of PRN's. We prove that, for any initial state, their pdf 
becomes Maxwellian in the N — > oo limit, after an infinite 
number of two-particle "collisions" take place. To this 
end, we first prove that our system is ergodic @Jl5| . The 
proof is not exceedingly difficult because our system is 
not deterministic. We also study its output as a function 
of N, and establish useful criteria for its implementation. 
Correlation test results are also reported. 

For the motivation, consider numbers v\, V2 ■ ■ ■ Vn, 
placed in N computer registers, analogous to velocities of 
N particles that make up a closed classical system in ID. 
Pairs of registers 1 and j, say, selected at random with- 
out bias, are to "interact" somehow, conserving quantity 
vf + v?. By analogy with the approach to equilibrium 
(i.e., to Maxwell's velocities distribution) that is believed 
to take place in Statistical Physics, we expect that suffi- 
cient number of iterations will lead to an approximately 
Gaussian pdf of register values, from which the desired 
PRN's may be drawn. (See also Ref. go|.) We define 
below the simplest interaction we can think of in order 
that (1) implementation on a computer be very fast, and 
(2) that we may be able to prove that a Gaussian pdf 
does indeed ensue. 

Before the algorithm is implemented, all TV registers 



must be initialized to, say, v l — 1 for all 1 satisfying 
1 < 1 < N, or all v t may be read from a set of N register 
values saved from a previous computer run, which we 
assume to fulfill = N. Let U(1,N), J7 8 (l,/V) be 

unbiased integer random variables, both in the interval 
[1,/V], except that U l cannot equal 1. The algorithm 
follows: 



i = U(l,N); j = Ui(l,N); 
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The updated value of v t , from Eq. (2), is used in Eq. (3). 
After an initial warm up phase (see below), v t and v 3 may 
be drawn each time transformation (1-3) is applied. They 
are two independent PRN's, each one with an approxi- 
mately Gaussian pdf, with (v % ) = and (vf) = 1 for all 1, 
if N is sufficiently large (see below). Transformation (1- 
3) may be thought of as a rotation of ±7r/4 with respect 
to a randomly chosen tj plane (+ and — signs are for the 
two possible index orderings, ij and ji). Thus, quantity 
^2 vf is conserved. Frequencies of events from sequences 
of 10 6 , 10 8 and 10 10 PRN's generated with transforma- 
tion (1-3), with N = 1024, are exhibited in Figj|. 
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FIG. 1. Number n{v) of PRN's generated within v — Aw/2 
and v + Av/2, for Aw = 0.1. Datapoins shown as •, □, and 
■ follow from sequences of 1.9 x 10 10 , 1.8 x 10 8 , and 2 x 10 6 
PRN's, respectively. 



We first explain why PRN's genered by transforma- 
tion (1-3) are expected to be normally distributed. Let 
P„(v) be the probability density at v = [v\, V2, • ■ ■ , vn), 
after transformation (1-3) has been applied n times, on 
the (N — l)-dimensional spherical surface Sn-i, of ra- 
dius VN given by N = J2e=i N v e- Let the single 
register pdf p(v) be the n — > 00 limit of p n {v), where 
Pn(vi) = J P n (v) dv2dv3 . . .oIvn- We show further below 
that P n (v) — ► constant over Sn-i as n — > 00. It then 
follows by integration that, 



p(v) oc 1 



N 



(JV-3)/2 
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Clearly, p(v) — > Cexp(— v 2 /2) in the N — > oo limit, 
which is the desired result. 

We prove below, in three stages, that P«(v) does in- 
deed become homogeneous over spherical surface Sn-i, 
if N > 3, in the n — ► oo limit. We first prove P n (v) 
P rl (u) as n — ► oo if v and u are related. [From here on, 
we say that points v and u are related if succesive trans- 
formations (1-3) of v can lead to u.] We then prove that 
the system's "orbit" covers iSat_i densely [that is, that 
any point v G 5/v-i can be brought arbitrarily close to 
any other point u G 5/v-i by applying transformations 
(1-3) to v a sufficient number of times]. Then, the de- 
sired result follows easily. It may help to place the signif- 
icance of the proof that follows into proper perspective 
to note that if in Eq. (1) j *— Ui[l,N] is replaced by 
j = t + 1 modN , the system becomes then non-ergodic, 
as can be easily checked numerically. 

To start the proof, let kernel K(v, v') be defined by 
P>+i(v) = J X(v,v')P„(v') dv', and let 



Fn = /{^+l(v)-P^(v)}dv. 



(0.5) 



Note first that F n < implies that P n+ i(v) is more uni- 
form than P„(v), in the sense that J dv [P n+1 (v) — P] 2 < 
/ civ [P n (v) — P] 2 , where P = 1/ j dv. It follows from 
the definition of if(v, v') that 



F n = J dv{[J dviif(v,vi)P n (vi 



)] 2 -P 2 (v)}. (0.6) 



Making use of the detailed balance condition, if(v, v') = 
K(v',v), which our sys tem satisfies, and the relation 
J dv K(v,v') — 1, Eq. (CK6) can be cast into, 



F n = ~fdvf dv x J dv 2 Q(v,vi, v a ), (0.7) 

where, Q = K(v,v 1 )K(v,v 2 )[P n (v 1 ) - P n (v 2 )} 2 . There- 
fore, in the n — > oo limit, P n (v) becomes constant over 
each set in Sn—i within which any two points v,u are 
related. 

We now prove that the system's orbit covers Sjy-i 
densely. Let Hn be the group of transformations in 
N dimensions defined by Eqs. (1-3). We first show 
that any rotation in 3D [that is, any element of 50(3)] 
can be approximated arbitrarily close by elements of H3. 
The proof is extended to higher dimensions by induc- 
tion. Note first that P3 does not belong to the set of 
finite rotation groups in 3D Jl6| , and is therefore an in- 
finite group. Let group SO (3) be covered by spheres of 
radius e/2 each. A finite number of them is sufficient, 
since the volume of 50(3) is finite fTij] . It follows that 
there must be at least one sphere with two elements of 
P3 in it, since P3 has an infinite number of elements. Let 
these two elements be r and s, and let 17(11, e) be element 
rs^ 1 of P3, which is a rotation by angle e about some 
undetermined u axis. We will build elements of P3 that 
are as near as desired to any given rotation. To this end, 



it is sufficient to show that it can be done for a set of 
infinitesimal generators of rotations |jisj| . One such set 
is made up of infinitesimal rotations about three linearly 
independent axes. Consider axes Ui, 112, and 113 that are 
obtained from u by rotations g(l,7r/2), <?(2,7r/2), and 
g(3, tt/2) about each one of the coordinate axes by angle 
7r/2. The correspondng infinitesimal rotations are given 
by [0, <7(u t ,e) = g(t,Tr/2)g(u,e)g- 1 (t,Tr/2). This con- 
eludes the proof for 3 dimensions. 

We now prove by induction that any element g(ij,a), 
for the rotation about plane ij, by angle a of the rotation 
group SO(N) can be approximated as nearly as desired 
by an element g of group Hn, for N > 3. By hypothesis, 
any g(ij, a), for i, j = 1, 2, . . . , N can be approximated 
by an element g of Hm- We show now that g(i N + 1, a), 
for 1 = 1, 2, . . . , TV, can also be approximated by elements 
of H^+i- We take g € Hn within distance e of g tJ (a). 
Now, since rotations preserve distances, it follows that 
g(t N + l,a) e SO{N + 1), given by g(t N + l,a) = 
g(i N + l, n/2)g(ij, a)<7 _1 (« N+l, tt/2) is within distance 
e of g' g H N+1 , given by g' = g(i N + /2)gg~ 1 (i N + 
1, 7r/2). This proves dense coverage in N > 3 dimensions. 
This is a stochastic generalization of Jacobi's theorem 
p"5| to more than two dimensions. 

To conclude the proof that P n (v) — > constant in the 
11 — > 00 limit, consider any two points V and U' as cen- 
ters of disks Dv an d 2\j', both of radius r, in Sn-i- 
Since the system's orbit covers Sn-i densely for N > 3, 
it follows that a point U that is related to V exists arbi- 
trarily close to U'. Consider now disks 2?v and T>\j. The 
fact that there exists at least one sequence of rotations in 
Hn that take V into U implies that there exists at least 
one single rotation g in Hn that transforms V into U. 
Since g is a rotation, it transforms 2?v rigidly into T>\j. 
It follows that J dv Pjy (v) over IV equals J du P/v (u) 
over T>u. Since r is arbitrary, and V and U' are any 
two points in iS/v-i, it follows that P(v) is constant over 
Sn-i (except, perhaps, on a set of measure zero). This 
is the desired result. Ergodicity follows [jl5f . 

We next address the following practical issues: (1) how 
good an approximation to a Gaussian pdf of PRN's is 
achieved with a necessarily finite set of N registers; (2) 
how long must the warm up phase be. 

It is convenient to rewrite Eq. (0.4) as follows, 



I \ -v 2 II q N (v)/N 

p[v) oc e ' e a ;/ . 



(0.8) 



where g N (v) = v 2 (3 - v 2 /2)/2 + 0(1/N). N^gnty) is 
approximately the fractional deviation, Sp(v)/p(v), from 
Gaussian form if Sp(v)/p(v) <C 1. We have checked this 
behavior numerically. Clearly, the number of registers 
N that must be used increases with the number M of 
PRN's one intends to generate. This is because the value 
of the largest PRN generated increases, on the average, 
with M. More precisely, the value of v beyond which 
PRN's are only generated with probability q is approxi- 
mate ly g iven by v 2 m 2hi(M/vq). Now, it follows from 
Eq. (3.8) that the fractional error SP/P in the probabil- 
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ity density at v is approximately N v (3 — v /2)/2 for 
very large N. (It is pointless to require this error to be 
too small since a PRN is expected to be generated be- 
yond x with a small probability q.) It then follows that 
[\n(M/qv)} 2 < N5P/P must be satisfied by N. Thus, 
approximately 10 4 registers are sufficient in order to gen- 
erate as many as 10 15 PRN's, with a roughly 10% error in 
the probability for the largest PRN in the sequence. For 
results obtained from a sequence of 10 10 PRN's generated 
with 1024 registers, see Fig. 1. 
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FIG. 2. Number n(v) of PRN's within v - Av/2 and 
v + Av/2, for Av — 0.1, starting from initial conditions v, = 1, 
for all i £ [1,N], after transformation (1-3) is iterated 2n p N 
times (that is, after each register interacts, on the average, 
2n p times). The □ and ♦ stand for n p = 2, 10, respectively, 
for TV = 1024. The O and •, and o stand for n v = 2, 4, and 10 
respectively, for N = 1048576. The two straight lines stand 
for C exp(— v 2 /2) for two values of C. 



in m succesively generated PRN's v±, V2, ■ ■ ■ v m , for m = 
3, 4, . . . , 6, performing a chi-square isotropy test over the 
corresponding m-dimensional space. An m-tuple v = 
vi,V2, ■ ■ ■ , v m was said to belong to the i—th cone, of 1024 
randomly oriented cones with axes Wi,W2, . . . , W1024, if 
0.99 < v.w, < 1. No significant deviations from isotropy 
were observed for 10 6 generated m-tuples. 
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FIG. 3. Average energy per spin, obtained from MC sim- 
ulations using Wolff's algorithm, versus the inverse of the 
number of registers used for the generation of PRN's with 
Gaussian pdf's. •, ■ and A stand for data points that follow 
from feeding our algorithm with the following uniform PRN 
generators: ggl, R(250,103,xor), and Ran3, respectively. Un- 
acceptable energy values that have been obtained in Refs. [3] 
using R(250,103,xor), and Ran3 are also shown as bars next 
to the y-axis. 



Our algorithm must be applied a number n p N of 
times before it is ready for use unless all v t are initial- 
ized with "equilibrium" values (stored from some previous 
computer run). The distribution of all register values 
then evolves towards equilibrium, as illustrated in Fig. 
2. Deviations from equilibrium are statistically insignif- 
icant for n p > 2 and N = 1024, and for n p > 4 and 
N = 1 048 576. Since n p is expected to increase as IniV, 
n p = 8 should provide ample warm up for any forseeable 
applications. 

The number of PRN's that must be generated be- 
fore each PRN in sequence v\,V2, ■ ■ ■ , vn returns within 
distance r from its initial value is exponential in N. 
More specifically, we estimate it to be (t / \f~N){\ / r) N 
for N ^S> 1, where r is the period of the algorithm used 
to select 1 and j in Eq. (1). The estimation is based on 
P n (v) — > constant over S^-i as n — > 00. Thus, an effec- 
tively infinite recurrence time follows for any reasonable 
value of N. 

Correlations between a finite number of PRN's clearly 
vanish as N — > 00, since 1 and j in Eq. (1) are supposedly 
independent PRN's. We have searched for correlations 



Implementation of Wolff's algorithm [[L9| in MC calcu- 
lations of the Ising model's critical behavior is a demand- 
ing test that some well known uniform PRN generators 
have failed S. Large clusters are then flipped as a whole, 
and this tests correlations in very long sequences. We 
have used normal PRN's generated by our algorithm as 
input into a MC simulation of an Ising system of 16 x 16 
spins at the critical temperature. [For that, we note that 
vj + v 2 > 2x as often as u > exp(— x) if v t and v 3 (it) 
are PRN's with Gaussian (uniform) pdf's, respectively] 
The energy obtained is shown in Fig. 3 as a function of 
the number of registers N. The following uniform PRN 
algorithms were used to select 1 and j in Eq. (1): ggl 
@, R(250, 103,xor) |||, and RAN3 @. We tried the 
latter two algorithms, which have been shown to lead by 
themselves to unacceptable results for the Ising model [[| , 
in order to test our algorithm's robustness. The results 
shown in Fig. 3 are gratifying. 

Similarly, the specific heat c and magnetization m fluc- 
tuations data points obtained follow approximately the 
relations c ~ c + 8A/N, and {{5m) 2 ) ~ xo + 33/7V, re- 
spectively, where Co = 1.497(1) and \o — 0.5454(2), in 
agreement with the known exact values p|,^l[. 
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Double precision is recommended. It prevents excesive 
drift of the sum vl away from its assigned value. Even 
then, single precision accuracy is to be expected at the 
end of a sequence of some 10 16 PRN's, unless the sum is 
normalized several times during the run. 

In summary, we have shown that implementation of 
Eqs. (1-3) provides a source of PRN's with an approx- 
imately Gaussian pdf. Some 10 4 registers (molecules) 
are sufficient for some purposes, but up to 10 5 or more 
may be necessary for more demanding tasks. (Having to 
make a decision about the number of registers to be used 
may sometimes be an unwelcomed task. On the other 
hand, it is a virtue of the algorithm, that one can con- 
trol, through the value of N, how close the output is to 
be from sequences of truly independent random numbers 
with Gaussian pdf's.) Initial warm ups for arbitrary ini- 
tial conditions are necessary; it is sufficient to let each 
register initially interact an average number of, say, 8 
times. The system's recurrence time was shown to be 
exponential in N, and therefore effectively infinite. Its 
behavior appears to be robust. The proposed algorithm 
runs an order of magnitude faster on computers than the 
most often used Box-Muller method ^,|). For a fortran 
code of our algorithm or other questions, please write 
JFF@Pipe.Unizar.Es. 
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